library(tidyverse)
library(purrrlyr)
library(lubridate)
library(themebg)
library(plotly)
library(broom)

x <- dirr::get_rds("../data/final")

se_line <- list(color = 'rgba(7, 164, 181, 0.05)')
se_fill <- 'rgba(7, 164, 181, 0.2)'

bolus <- hep_bolus_initiation %>%
    distinct(millennium.id) %>%
    mutate(bolus = TRUE)
df <- ptt_run %>%
    left_join(bolus, by = "millennium.id") %>%
    left_join(hep_protocol[c("millennium.id", "order")], by = "millennium.id") %>%
    dmap_at("bolus", ~coalesce(.x, FALSE)) %>%
    filter(!is.na(duration_heparin),
           duration_hypothermia > -24,
           duration_hypothermia < 72,
           order %in% c("ACS", "DVT", "Stroke")) %>%
    left_join(hep_drip[c("millennium.id", "start.datetime", "stop.datetime")], by = "millennium.id") %>%
    filter(lab.datetime >= start.datetime - hours(4),
           lab.datetime <= stop.datetime + hours(4))
hep_run %>%
    filter(drip.count == 1,
           order %in% c("ACS", "DVT", "Stroke")) %>%
    ggplot(aes(x = run.time, y = med.rate, color = order)) +
    geom_point(shape = 1) +
    # geom_line() +
    geom_smooth() +
    xlim(c(0, 48)) +
    # facet_wrap(~ order, ncol = 1) +
    theme_bg()
hep_run %>%
    filter(drip.count == 1,
           order %in% c("ACS", "DVT", "Stroke")) %>%
    ggplot(aes(x = duration_hypothermia, y = med.rate, color = order)) +
    geom_point(shape = 1) +
    # geom_line() +
    geom_smooth() +
    xlim(c(-24, 72)) +
    theme_bg()
ptt_run %>%
    left_join(bolus, by = "millennium.id") %>%
    left_join(hep_protocol[c("millennium.id", "order")], by = "millennium.id") %>%
    dmap_at("bolus", ~coalesce(.x, FALSE)) %>%
    filter(order %in% c("ACS", "DVT")) %>%
    ggplot(aes(x = duration_heparin, y = lab.result, color = order)) +
    geom_point(shape = 1) +
    geom_smooth() +
    xlim(c(-12, 96)) +
    theme_bg()
ptt_run %>%
    # mutate(on_heparin = lab.datetime >= heparin.start) %>%
    left_join(bolus, by = "millennium.id") %>%
    dmap_at(c("bolus", "on_heparin"), ~coalesce(.x, FALSE)) %>%
    filter(!is.na(duration_heparin)) %>%
    ggplot(aes(x = duration_hypothermia, y = lab.result, color = bolus)) +
    # geom_point(aes(shape = on_heparin)) +
    geom_vline(xintercept = 0, color = "light grey") +
    geom_point(shape = 1) +
    geom_smooth() +
    scale_shape_manual(values = c(1, 2, 3)) +
    # xlim(c(-12, 96)) +
    scale_x_continuous(breaks = seq(-24, 120, 12), limits = c(-12, 96)) +
    theme_bg()
df <- ptt_run %>%
    left_join(bolus, by = "millennium.id") %>%
    left_join(hep_protocol[c("millennium.id", "order")], by = "millennium.id") %>%
    dmap_at("bolus", ~coalesce(.x, FALSE)) %>%
    filter(!is.na(duration_heparin),
           duration_hypothermia > -24,
           duration_hypothermia < 72,
           order %in% c("ACS", "DVT", "Stroke")) %>%
    left_join(hep_drip[c("millennium.id", "start.datetime", "stop.datetime")], by = "millennium.id") %>%
    filter(lab.datetime >= start.datetime - hours(4),
           lab.datetime <= stop.datetime + hours(4))

d_acs <- filter(df, order == "ACS")
d_dvt <- filter(df, order == "DVT")
d_cva <- filter(df, order == "Stroke")

m_acs <- loess(lab.result ~ duration_hypothermia, data = df, subset = order == "ACS")
m_dvt <- loess(lab.result ~ duration_hypothermia, data = df, subset = order == "DVT")
m_cva <- loess(lab.result ~ duration_hypothermia, data = df, subset = order == "Stroke")

plot_ly(x = ~duration_hypothermia) %>%
    add_markers(y = ~lab.result, data = df, split = ~order, marker = list(symbol = "circle-open")) %>%
    # add_markers(y = ~lab.result, data = d_acs, marker = list(symbol = "circle-open")) %>%
    # add_markers(y = ~lab.result, data = d_dvt, marker = list(symbol = "circle-open")) %>%
    add_lines(y = ~fitted(m_acs), data = d_acs, name = "ACS", legendgroup = "acs") %>%
    add_lines(y = ~fitted(m_dvt), data = d_dvt, name = "DVT", legendgroup = "dvt") %>%
    add_lines(y = ~fitted(m_cva), data = d_cva, name = "Stroke", legendgroup = "cva") %>%
    add_ribbons(data = augment(m_acs),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = list(color = 'rgba(7, 164, 181, 0.05)'),
                fillcolor = 'rgba(7, 164, 181, 0.2)',
                showlegend = FALSE,
                legendgroup = "acs") %>%
    add_ribbons(data = augment(m_dvt),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = list(color = 'rgba(7, 164, 181, 0.05)'),
                fillcolor = 'rgba(7, 164, 181, 0.2)',
                showlegend = FALSE, 
                legendgroup = "dvt") %>%
        add_ribbons(data = augment(m_cva),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = list(color = 'rgba(7, 164, 181, 0.05)'),
                fillcolor = 'rgba(7, 164, 181, 0.2)',
                showlegend = FALSE, 
                legendgroup = "cva")

PTT change over time from hypothermia initiation by heparin protocol

# hep_bolus_initial %>%
#     left_join(hep_start, by = "millennium.id") %>%
#     left_join(hypothermia_start, by = "millennium.id") %>%
#     filter(!is.na(hypothermia_start)) %>%
#     mutate(time_heparin = as.numeric(difftime(med.datetime, heparin.start, units = "hours")),
#            time_hypothermia = as.numeric(difftime(med.datetime, hypothermia_start, units = "hours"))) %>%
#     filter(time_hypothermia >= -24, time_hypothermia <= 72) %>%
#     ggplot(aes(x = time_hypothermia, y = med.dose)) +
#     geom_point(shape = 1) +
#     theme_bg()

Univariate Comparisons

d <- data_wt_avg %>%
    left_join(patients, by = "millennium.id") %>%
    filter(!is.na(temp.time.wt.avg),
           !is.na(ptt.time.wt.avg))

m <- lm(ptt.time.wt.avg ~ temp.time.wt.avg, data = d)

plot_ly(d, x = ~temp.time.wt.avg) %>%
    add_markers(y = ~ptt.time.wt.avg, split = ~group, marker = list(symbol = "circle-open")) %>%
    add_lines(y = ~fitted(m), showlegend = FALSE) %>%
    add_ribbons(data = augment(m),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = list(color = 'rgba(7, 164, 181, 0.05)'),
                fillcolor = 'rgba(7, 164, 181, 0.2)',
                showlegend = FALSE)

Time-Weighted Average of Temperature vs. PTT

mod <- lm(ptt.time.wt.avg ~ temp.time.wt.avg, data_wt_avg) 
    
mod %>%
    glance() %>%
    knitr::kable(caption = "Linear Regression Model for Time-Weighted Average PTT by Temperature")
Linear Regression Model for Time-Weighted Average PTT by Temperature
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
0.0026356 -0.0029053 25.27819 0.4756583 0.4912837 2 -845.0907 1696.181 1705.793 115017.6 180
mod %>%
    tidy() %>%
    knitr::kable(caption = "Coefficients for PTT by Temperature")
Coefficients for PTT by Temperature
term estimate std.error statistic p.value
(Intercept) 155.3361788 134.696583 1.1532303 0.2503450
temp.time.wt.avg -0.9894958 1.434718 -0.6896798 0.4912837
d <- temp_ptt_full %>%
    left_join(patients, by = "millennium.id") %>%
    filter(!is.na(vital.result),
           !is.na(lab.result),
           vital.result > 85)

d_study <- filter(d, group == "Study")
d_contr <- filter(d, group == "Control")

# m <- lm(lab.result ~ vital.result, data = d)
m_study <- lm(lab.result ~ vital.result, data = d, subset = group == "Study")
m_contr <- lm(lab.result ~ vital.result, data = d, subset = group == "Control")

plot_ly(x = ~vital.result) %>%
    add_markers(data = d, 
                y = ~lab.result, 
                split = ~group, 
                marker = list(symbol = "circle-open")) %>%
    add_lines(y = ~fitted(m_study), 
              data = d_study, 
              legendgroup = "study", 
              name = "Study") %>%
    add_lines(y = ~fitted(m_contr), 
              data = d_contr, 
              legendgroup = "control", 
              name = "Control") %>%
    add_ribbons(data = augment(m_study),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = se_line,
                fillcolor = se_fill,
                legendgroup = "study",
                showlegend = FALSE) %>%
    add_ribbons(data = augment(m_contr),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = se_line,
                fillcolor = se_fill,
                legendgroup = "control",
                showlegend = FALSE)

Temperature vs. PTT

mod <- lm(lab.result ~ vital.result, temp_ptt_full)

mod %>%
    glance() %>%
    knitr::kable(caption = "Linear Regression Model for PTT by Temperature")
Linear Regression Model for PTT by Temperature
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
0.0385579 0.0383575 44.42905 192.4199 0 2 -25020.59 50047.19 50066.62 9470965 4798
mod %>%
    tidy() %>%
    knitr::kable(caption = "Coefficients for PTT by Temperature")
Coefficients for PTT by Temperature
term estimate std.error statistic p.value
(Intercept) 335.309830 19.0040949 17.64408 0
vital.result -2.780202 0.2004247 -13.87155 0
d <- temp_ptt_full %>%
    left_join(patients, by = "millennium.id") %>%
    filter(!is.na(vital.result),
           !is.na(med.rate),
           vital.result > 85)

d_study <- filter(d, group == "Study")
d_contr <- filter(d, group == "Control")

m_study <- lm(med.rate ~ vital.result, data = d, subset = group == "Study")
m_contr <- lm(med.rate ~ vital.result, data = d, subset = group == "Control")

plot_ly(x = ~vital.result) %>%
    add_markers(data = d, 
                y = ~med.rate, 
                split = ~group, 
                marker = list(symbol = "circle-open")) %>%
    add_lines(y = ~fitted(m_study), 
              data = d_study, 
              legendgroup = "study", 
              name = "Study") %>%
    add_lines(y = ~fitted(m_contr), 
              data = d_contr, 
              legendgroup = "control", 
              name = "Control") %>%
    add_ribbons(data = augment(m_study),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = se_line,
                fillcolor = se_fill,
                legendgroup = "study",
                showlegend = FALSE) %>%
    add_ribbons(data = augment(m_contr),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = se_line,
                fillcolor = se_fill,
                legendgroup = "control",
                showlegend = FALSE)

Temperature vs. Heparin

temp_ptt_full %>%
    ungroup() %>%
    # select(-study_id) %>%
    select_if(is.numeric) %>%
    cor(use = "pairwise.complete.obs") %>%
    tidy() %>%
    knitr::kable(caption = "Correlation Analysis")
Correlation Analysis
.rownames vital.result lab.result med.rate
vital.result 1.0000000 -0.1963616 0.2709087
lab.result -0.1963616 1.0000000 -0.1404534
med.rate 0.2709087 -0.1404534 1.0000000
d <- data_wt_avg %>%
    left_join(patients, by = "millennium.id") %>%
    filter(!is.na(hep.time.wt.avg),
           !is.na(ptt.time.wt.avg))

m <- lm(ptt.time.wt.avg ~ hep.time.wt.avg, data = d)

plot_ly(d, x = ~hep.time.wt.avg) %>%
    add_markers(y = ~ptt.time.wt.avg, split = ~group, marker = list(symbol = "circle-open")) %>%
    add_lines(y = ~fitted(m), showlegend = FALSE) %>%
    add_ribbons(data = augment(m),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = list(color = 'rgba(7, 164, 181, 0.05)'),
                fillcolor = 'rgba(7, 164, 181, 0.2)',
                showlegend = FALSE)

Time-Weighted Average of Heparin Dose vs. PTT

mod <- lm(ptt.time.wt.avg ~ hep.time.wt.avg, data_wt_avg) 
    
mod %>%
    glance() %>%
    knitr::kable(caption = "Linear Regression Model for Time-Weighted Average PTT by Heparin Dose")
Linear Regression Model for Time-Weighted Average PTT by Heparin Dose
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
0.1084269 0.1007409 24.07506 14.10711 0.0002717 2 -541.805 1089.61 1097.922 67234.59 116
mod %>%
    tidy() %>%
    knitr::kable(caption = "Coefficients for PTT by Heparin Dose")
Coefficients for PTT by Heparin Dose
term estimate std.error statistic p.value
(Intercept) 87.807266 5.1864449 16.930145 0.0000000
hep.time.wt.avg -1.922917 0.5119666 -3.755943 0.0002717
d <- temp_ptt %>%
    left_join(patients, by = "millennium.id") %>%
    filter(!is.na(ptt),
           !is.na(hep),
           temp > 85)

# m <- loess(ptt ~ hep, data = d)
m <- lm(ptt ~ hep, data = d)

plot_ly(d, x = ~hep) %>%
    add_markers(y = ~ptt, split = ~group, marker = list(symbol = "circle-open")) %>%
    add_lines(y = ~fitted(m), showlegend = FALSE) %>%
    add_ribbons(data = augment(m),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = list(color = 'rgba(7, 164, 181, 0.05)'),
                fillcolor = 'rgba(7, 164, 181, 0.2)',
                showlegend = FALSE)

Heparin Dose vs. PTT

mod <- lm(ptt ~ hep, temp_ptt)

mod %>%
    glance() %>%
    knitr::kable(caption = "Linear Regression Model for PTT by Heparin Dose")
Linear Regression Model for PTT by Heparin Dose
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
0.0373575 0.0343303 47.44639 12.34071 0.0005074 2 -1688.129 3382.259 3393.564 715869 318
mod %>%
    tidy() %>%
    knitr::kable(caption = "Coefficients for PTT by Heparin Dose")
Coefficients for PTT by Heparin Dose
term estimate std.error statistic p.value
(Intercept) 63.741656 6.2418938 10.211910 0.0000000
hep 2.191538 0.6238482 3.512935 0.0005074
d <- data_wt_avg %>%
    left_join(patients, by = "millennium.id") %>%
    filter(!is.na(hep.time.wt.avg),
           !is.na(temp.time.wt.avg))

m <- lm(hep.time.wt.avg ~ temp.time.wt.avg, data = d)

plot_ly(d, x = ~temp.time.wt.avg) %>%
    add_markers(y = ~hep.time.wt.avg, split = ~group, marker = list(symbol = "circle-open")) %>%
    add_lines(y = ~fitted(m), showlegend = FALSE) %>%
    add_ribbons(data = augment(m),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = list(color = 'rgba(7, 164, 181, 0.05)'),
                fillcolor = 'rgba(7, 164, 181, 0.2)',
                showlegend = FALSE)

Time-Weighted Average of Temperature vs. Heparin Dose

mod <- lm(hep.time.wt.avg ~ temp.time.wt.avg, data_wt_avg) 
    
mod %>%
    glance() %>%
    knitr::kable(caption = "Linear Regression Model for Time-Weighted Average Heparin Dose by Temperature")
Linear Regression Model for Time-Weighted Average Heparin Dose by Temperature
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
0.0480119 0.0398051 4.260029 5.850261 0.0171294 2 -337.4407 680.8815 689.1935 2105.15 116
mod %>%
    tidy() %>%
    knitr::kable(caption = "Coefficients for Heparin Dose by Temperature")
Coefficients for Heparin Dose by Temperature
term estimate std.error statistic p.value
(Intercept) -55.2210811 26.6201481 -2.074409 0.0402537
temp.time.wt.avg 0.6854836 0.2834063 2.418731 0.0171294
d <- temp_ptt %>%
    left_join(patients, by = "millennium.id") %>%
    filter(!is.na(temp),
           !is.na(hep),
           temp > 85)

m <- lm(hep ~ temp, data = d)

plot_ly(d, x = ~temp) %>%
    add_markers(y = ~hep, split = ~group, marker = list(symbol = "circle-open")) %>%
    add_lines(y = ~fitted(m), showlegend = FALSE) %>%
    add_ribbons(data = augment(m),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = list(color = 'rgba(7, 164, 181, 0.05)'),
                fillcolor = 'rgba(7, 164, 181, 0.2)',
                showlegend = FALSE)

Temperature vs. Heparin Dose

mod <- lm(hep ~ temp, temp_ptt)

mod %>%
    glance() %>%
    knitr::kable(caption = "Linear Regression Model for Heparin Dose by Temperature")
Linear Regression Model for Heparin Dose by Temperature
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
0.0077065 0.0045861 4.248453 2.469699 0.1170546 2 -915.9547 1837.909 1849.214 5739.693 318
mod %>%
    tidy() %>%
    knitr::kable(caption = "Coefficients for Heparin Dose by Temperature")
Coefficients for Heparin Dose by Temperature
term estimate std.error statistic p.value
(Intercept) 21.5997124 7.9846045 2.705170 0.0071949
temp -0.1337619 0.0851158 -1.571528 0.1170546

Multiple Variable Models

mod <- lm(ptt.time.wt.avg ~ temp.time.wt.avg + hep.time.wt.avg, data_wt_avg) 
    
mod %>%
    glance() %>%
    knitr::kable(caption = "Linear Regression Model for Time-Weighted Average PTT by Temperature and Heparin Dose")
Linear Regression Model for Time-Weighted Average PTT by Temperature and Heparin Dose
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
0.1109787 0.0955174 24.14488 7.177862 0.0011547 3 -541.6359 1091.272 1102.355 67042.16 115
mod %>%
    tidy() %>%
    knitr::kable(caption = "Coefficients for PTT by Temperature and Heparin Dose")
Coefficients for PTT by Temperature and Heparin Dose
term estimate std.error statistic p.value
(Intercept) 176.0334300 153.6499874 1.1456781 0.2543069
temp.time.wt.avg -0.9458452 1.6462895 -0.5745315 0.5667304
hep.time.wt.avg -1.8566694 0.5262396 -3.5281826 0.0006024
d <- temp_ptt_full %>%
    left_join(hep_bolus_initiation[c("millennium.id", "med.dose")], by = "millennium.id") %>%
    mutate(bolus = !is.na(med.dose))

mod <- lm(lab.result ~ vital.result + med.rate + order + bolus, d)

mod %>%
    glance() %>%
    knitr::kable(caption = "Linear Regression Model for PTT by Temperature and Heparin Dose")
Linear Regression Model for PTT by Temperature and Heparin Dose
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
0.1720592 0.1704976 43.42129 110.177 0 7 -16541.86 33099.72 33148.26 5997486 3181
mod %>%
    tidy() %>%
    knitr::kable(caption = "Coefficients for PTT by Temperature and Heparin Dose")
Coefficients for PTT by Temperature and Heparin Dose
term estimate std.error statistic p.value
(Intercept) 548.6496384 23.9237290 22.9332826 0.0000000
vital.result -4.9427775 0.2544804 -19.4230194 0.0000000
med.rate -0.3153441 0.1210017 -2.6061132 0.0092003
orderDVT 27.0595968 2.0263971 13.3535509 0.0000000
orderFixed 2.0114091 7.7920944 0.2581346 0.7963197
orderStroke -9.4901118 3.1242578 -3.0375572 0.0024044
bolusTRUE 9.3355311 1.7299284 5.3964842 0.0000001

Chelsea’s Data

plot_ly(ck_heparin, x = ~as.numeric(duration), y = ~Value, split = ~Control) %>%
    add_markers()
plot_ly(ck_ptt, x = ~duration_ptt, y = ~Value, split = ~Control) %>%
    add_markers(marker = list(symbol = "circle-open"))
d <- ck_data %>%
    filter(!is.na(Temperature),
           !is.na(heparin))
           # heparin > 0,
           # heparin < 4000)

m <- lm(heparin ~ Temperature, data = d)

plot_ly(d, x = ~Temperature) %>%
    add_markers(y = ~heparin, split = ~Control, marker = list(symbol = "circle-open")) %>%
    add_lines(y = ~fitted(m), showlegend = FALSE) %>%
    add_ribbons(data = augment(m),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = list(color = 'rgba(7, 164, 181, 0.05)'),
                fillcolor = 'rgba(7, 164, 181, 0.2)',
                showlegend = FALSE)

Temperature vs. Heparin Dose

d <- ck_data %>%
    filter(!is.na(Temperature),
           !is.na(PTT))

d_study <- filter(d, Control == "Study")
d_ctrl <- filter(d, Control == "Control")

m_study <- lm(PTT ~ Temperature, data = d_study)
m_ctrl <- lm(PTT ~ Temperature, data = d_ctrl)

plot_ly(x = ~Temperature) %>%
    add_markers(data = d, y = ~PTT, split = ~Control, marker = list(symbol = "circle-open")) %>%
    add_lines(data = d_study, y = ~fitted(m_study), name = "Study", legendgroup = "study") %>%
    add_lines(data = d_ctrl, y = ~fitted(m_ctrl), name = "Control", legendgroup = "control") %>%
    add_ribbons(data = augment(m_study),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = list(color = 'rgba(7, 164, 181, 0.05)'),
                fillcolor = 'rgba(7, 164, 181, 0.2)',
                legendgroup = "study",
                showlegend = FALSE) %>%
    add_ribbons(data = augment(m_ctrl),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = list(color = 'rgba(7, 164, 181, 0.05)'),
                fillcolor = 'rgba(7, 164, 181, 0.2)',
                legendgroup = "control",
                showlegend = FALSE) 

Temperature vs. PTT